###########################################################################################################################################################
# Replication Lewandowsky/Jankowski: Sympathy for the devil.
# Forthcoming: EPSR
# Contact: michael.jankowski@uol.de
###########################################################################################################################################################

library(tidyverse)
library(Zelig)

illiberal <- readRDS("finaldata.rds")

##################################################
# Logit models
##################################################

logit <- zelig(selected ~ closeness_advantage_illiberal_w +
                 populism +
                 satdem +
                 educ +
                 auth +
                 pluralism +
                 leftright +
                 trust_reg +
                 polint_cat +
                 gender +
                 age +
                 party,
               data = illiberal,
               model = "logit")

logit2 <- zelig(selected ~ closeness_advantage_illiberal_w*populism +
                  satdem +
                  educ +
                  auth +
                  pluralism +
                  leftright +
                  trust_reg +
                  polint_cat +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")

logit3 <- zelig(selected ~ closeness_advantage_illiberal_w*auth+
                  satdem +
                  populism +
                  educ +
                  pluralism +
                  leftright +
                  trust_reg +
                  polint_cat +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")

logit4 <- zelig(selected ~ closeness_advantage_illiberal_w*auth*populism +
                  satdem +
                  educ +
                  pluralism +
                  leftright +
                  trust_reg +
                  polint_cat +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")


logit5 <- zelig(selected ~ closeness_advantage_illiberal_w*leftright +
                  auth + populism +
                  satdem +
                  educ +
                  pluralism +
                  trust_reg +
                  polint_cat +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")


logit6 <- zelig(selected ~ closeness_advantage_illiberal_w*leftright*I(leftright^2)  +
                  auth + populism +
                  satdem +
                  educ +
                  pluralism +
                  trust_reg +
                  polint_cat +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")

logit7 <- zelig(selected ~ closeness_advantage_illiberal_w*polint_cat + 
                  leftright +
                  satdem +
                  auth + 
                  populism +
                  educ +
                  pluralism +
                  trust_reg +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")

logit8 <- zelig(selected ~ closeness_advantage_illiberal_w*satdem +
                  polint_cat + 
                  leftright +
                  auth + 
                  populism +
                  educ +
                  pluralism +
                  trust_reg +
                  gender +
                  age +
                  party,
                data = illiberal,
                model = "logit")

logit9 <- zelig(selected ~ closeness_advantage_illiberal_w*party +
                  satdem +
                  polint_cat + 
                  leftright +
                  auth + 
                  populism +
                  educ +
                  pluralism +
                  trust_reg +
                  gender +
                  age,
                data = illiberal,
                model = "logit")

illiberal$party_grob <- as.factor(illiberal$party_grob)

logit10 <- zelig(selected ~ closeness_advantage_illiberal_w*party_grob +
                  satdem +
                  polint_cat + 
                  leftright +
                  auth + 
                  populism +
                  educ +
                  pluralism +
                  trust_reg +
                  gender +
                  age,
                data = illiberal,
                model = "logit")

##################################################
# Visulaize Model 1
##################################################

range_closeness <- seq(-2,2, length.out = 100)
range_breaks <- seq(-2,2,by = .5)

Model.close <-setx(logit, 
                   closeness_advantage_illiberal_w = range_closeness)

qoi <- sim(logit, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w) %>%
  summarise(expected_value = median(expected_value))

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = group),
            color = "purple", alpha =.05) +
  geom_line(data = median_val, color = "white") +
  theme_bw(base_size = 14) +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  lims(y = c(0,1))

ggsave(file = "Figure1.pdf", width = 7, height = 4.5)

##################################################
# Visulaize Model 2 (Populism)
##################################################

Model.close <- setx(logit2, 
                   closeness_advantage_illiberal_w = range_closeness,
                   populism = c(1,5))

qoi <- sim(logit2, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, populism) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, populism) %>%
  summarise(expected_value = median(expected_value))

preds$Populism <- ifelse(preds$populism == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))

median_val$Populism <- ifelse(median_val$populism == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Populism), 
                color = Populism),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Populism,
                                   color = Populism)) +
  geom_line(data = median_val, aes(group = Populism),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange")) +
  lims(y = c(0,1)) 

ggsave(file = "Figure2.pdf", width = 7, height = 4.5)

# First Difference

fd_fun <- function(range_closeness = 0){

Model.closea1 <- setx(logit2, 
                    closeness_advantage_illiberal_w = range_closeness,
                    populism = c(2))

Model.closea2 <- setx1(logit2, 
                     closeness_advantage_illiberal_w = range_closeness,
                     populism = c(5))

qoi <- Zelig::sim(logit2, 
           x = Model.closea1,
           x1 = Model.closea2,
           num = 1000)

fd.med <- median(qoi$sim.out$x1$fd[[1]])
fd.lo <- quantile(qoi$sim.out$x1$fd[[1]], prob=0.025)
fd.up <- quantile(qoi$sim.out$x1$fd[[1]], prob=0.975)
cat("fd.med = ", fd.med, "\n\n"); flush.console()

data.frame(med = fd.med, 
           lower = fd.lo, 
           upper = fd.up, 
           range_closeness = range_closeness)
}

fds <- map_df(seq(-2,2,by = .01), fd_fun)

ggplot(fds, aes(x = range_closeness,
                y = med,
                ymin = lower,
                ymax = upper)) + 
  geom_ribbon(alpha = .5) +
  geom_line() +
  theme_bw()
  
##################################################
# Visulaize Model 3 (Auth)
##################################################

Model.close <- setx(logit3, 
                    closeness_advantage_illiberal_w = range_closeness,
                    auth = c(1,5))

qoi <- sim(logit3, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, auth) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, auth) %>%
  summarise(expected_value = median(expected_value))

preds$Authoritarianism <- ifelse(preds$auth == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))

median_val$Authoritarianism <- ifelse(median_val$auth == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Authoritarianism), 
                color = Authoritarianism),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Authoritarianism,
                                   color = Authoritarianism)) +
  geom_line(data = median_val, aes(group = Authoritarianism),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange")) +
  lims(y = c(0,1))

ggsave(file = "Figure3.pdf", width = 7, height = 4.5)

##################################################
# Visualize Model 4 (Auth X Populsim)
##################################################

Model.close1 <- setx(logit4, 
                     closeness_advantage_illiberal_w = range_closeness,
                     auth = c(1,5),
                     populism = c(1))

Model.close2 <- setx(logit4, 
                     closeness_advantage_illiberal_w = range_closeness,
                     auth = c(1,5),
                     populism = c(5))

qoi <- sim(logit4, 
           x = Model.close1, 
           x1 = Model.close2)

preds <- zelig_qi_to_df(qoi)
preds <- preds %>% group_by(closeness_advantage_illiberal_w,
                            auth,
                            populism) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w,
           auth, 
           populism) %>%
  summarise(expected_value = median(expected_value))


preds$Authoritarianism <- ifelse(preds$auth == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))
median_val$Authoritarianism <- ifelse(median_val$auth == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Low", "High"))

preds$Populism <- ifelse(preds$populism == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Populism: Low", "Populism: High"))
median_val$Populism <- ifelse(median_val$populism == 1, 1, 2) %>%
  factor(levels = 1:2, labels = c("Populism: Low", "Populism: High"))

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Authoritarianism), 
                color = Authoritarianism),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Authoritarianism,
                                   color = Authoritarianism)) +
  geom_line(data = median_val, aes(group = Authoritarianism),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  facet_wrap(~ Populism) +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange")) +
  lims(y = c(0,1))

ggsave(file = "Figure4.pdf", width = 10, height = 6)

##################################################
# Visulaize Model 5 (LR)
##################################################

Model.close <- setx(logit5, 
                    closeness_advantage_illiberal_w = range_closeness,
                    leftright = c(1,5,11))

qoi <- sim(logit5, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, leftright) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, leftright) %>%
  summarise(expected_value = median(expected_value))

preds$Leftright <- ifelse(preds$leftright == 1, "Left",
                          ifelse(preds$leftright == 5, "Center", "Right")) %>%
  factor(labels = c("Left", "Center", "Right"))

median_val$Leftright <- ifelse(median_val$leftright == 1, "Left",
                          ifelse(median_val$leftright == 5, "Center", "Right")) %>%
  factor(labels = c("Left", "Center", "Right"))


ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Leftright), 
                color = Leftright),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Leftright,
                                   color = Leftright)) +
  geom_line(data = median_val, aes(group = Leftright),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange", "steelblue")) +
  lims(y = c(0,1)) +
  facet_wrap(~ Leftright) +
  guides(color = F)

ggsave(file = "Figure5.pdf", width = 10, height = 4.5)

##################################################
# Visulaize Model 6 (LR2)
##################################################

Model.close <- setx(logit6, 
                    closeness_advantage_illiberal_w = range_closeness,
                    leftright = c(1,5,11))

qoi <- sim(logit6, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, leftright) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, leftright) %>%
  summarise(expected_value = median(expected_value))

preds$Leftright <- ifelse(preds$leftright == 1, "Left",
                          ifelse(preds$leftright == 5, "Center", "Right")) %>%
  factor(labels = c("Left", "Center", "Right"))

median_val$Leftright <- ifelse(median_val$leftright == 1, "Left",
                               ifelse(median_val$leftright == 5, "Center", "Right")) %>%
  factor(labels = c("Left", "Center", "Right"))


ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Leftright), 
                color = Leftright),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Leftright,
                                   color = Leftright)) +
  geom_line(data = median_val, aes(group = Leftright),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange", "steelblue")) +
  lims(y = c(0,1)) +
  facet_wrap(~ Leftright) +
  guides(color = F)

ggsave(file = "Figure6.pdf", width = 10, height = 4.5)

##################################################
# Figure 8 (POl int) based on logit7
##################################################

Model.close <- setx(logit7, 
                    closeness_advantage_illiberal_w = range_closeness,
                    polint_cat = c("Low", "Moderate", "High"))

qoi <- sim(logit7, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, polint_cat) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, polint_cat) %>%
  summarise(expected_value = median(expected_value))

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, polint_cat), 
                color = polint_cat),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = polint_cat,
                                   color = polint_cat)) +
  geom_line(data = median_val, aes(group = polint_cat),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual(values = c("purple", "orange", "steelblue")) +
  lims(y = c(0,1)) +
  facet_wrap(~ polint_cat) +
  guides(color = F)

ggsave(file = "Figure8.pdf", width = 10, height = 4.5)

##################################################
# Visulaize Model 8 SATDEM --> Figure 9
##################################################

Model.close <- setx(logit8, 
                    closeness_advantage_illiberal_w = range_closeness,
                    satdem = c(1,5))

qoi <- sim(logit8, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, satdem) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, satdem) %>%
  summarise(expected_value = median(expected_value))

preds$Satdem <- ifelse(preds$satdem == 1, 
                       "Satisfied",
                       "Not Satisfied") %>%
  factor(labels = c("Satisfied", "Not Satisfied"))

median_val$Satdem <- ifelse(median_val$satdem == 1, 
                       "Satisfied",
                       "Not Satisfied") %>%
  factor(labels = c("Satisfied", "Not Satisfied"))


ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, Satdem), 
                color = Satdem),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = Satdem,
                                   color = Satdem)) +
  geom_line(data = median_val, aes(group = Satdem),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual("Sat. with Democracy in Germany", values = c("purple", "orange", "steelblue")) +
  lims(y = c(0,1))

ggsave(file = "Figure9.pdf", width = 10, height = 6)

##################################################
# Visulaize Model 9 Party Raw
##################################################

Model.close <- setx(logit9, 
                    closeness_advantage_illiberal_w = range_closeness,
                    party = sort(paste(unique(illiberal$party))))

qoi <- Zelig::sim(logit9, 
           x = Model.close,
           num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, party) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, party) %>%
  summarise(expected_value = median(expected_value))

median_val$party <- paste(median_val$party)
preds$party <- paste(preds$party)

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, party), 
                color = party),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = party,
                                   color = party)) +
  geom_line(data = median_val, aes(group = party),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual("Vote Choice", values = c("steelblue", 
                                               "darkgreen", 
                                               "blue",
                                               "purple", 
                                               "orange", 
                                               "darkred")) +
  lims(y = c(0,1)) +
  facet_wrap(~ party) -> partyplot

partyplot

ggsave(partyplot, file = "FigureA1.pdf", width = 10, height = 6)

##################################################
# Visulaize Model 10 PARYTGROB --> Figure 7
##################################################

Model.close <- setx(logit10, 
                    closeness_advantage_illiberal_w = range_closeness,
                    party_grob = sort(paste(unique(illiberal$party_grob))))

qoi <- Zelig::sim(logit10, 
                  x = Model.close,
                  num = 1000)

preds <- zelig_qi_to_df(qoi) %>%
  group_by(closeness_advantage_illiberal_w, party_grob) %>%
  mutate(group = 1:n())

median_val <- preds %>% 
  group_by(closeness_advantage_illiberal_w, party_grob) %>%
  summarise(expected_value = median(expected_value))

median_val$party_grob <- paste(median_val$party_grob)
preds$party_grob <- paste(preds$party_grob)

ggplot(data = preds, aes(x = closeness_advantage_illiberal_w,
                         y = expected_value)) +
  geom_line(aes(group = paste(group, party_grob), 
                color = party_grob),
            show.legend = F, alpha = .05) +
  geom_line(data = median_val, aes(group = party_grob,
                                   color = party_grob)) +
  geom_line(data = median_val, aes(group = party_grob),
            color = "white",
            show.legend = F) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Pr(Vote for Illiberal Candidate)") +
  xlab("Relative Weighted Policy Advantage") +
  scale_x_continuous(breaks = range_breaks) +
  scale_color_manual("Vote Choice", values = c("grey50", 
                                               "purple", 
                                               "blue")) +
  lims(y = c(0,1)) +
  facet_wrap(~ party_grob) -> party_grob_plot

party_grob_plot

ggsave(party_grob_plot, file = "Figure7.pdf", 
       width = 10, height = 4.5)

#####################################################
# Logit Table A1
#####################################################

texreg::texreg(list(logit,logit2,logit3,logit4,logit5,logit6,logit7,logit8,logit10),
               booktabs = TRUE)

#####################################################
# Table A2
#####################################################

library(psych)

model_data <- select(illiberal, 
                     closeness_advantage_illiberal_w,
                     populism,
                     satdem,
                     auth,
                     pluralism,
                     leftright,
                     trust_reg,
                     age,
                     gender,
                     educ,
                     polint_cat,
                     party) %>% na.omit()

describe(model_data) -> sum_cont
sum_cont$varname <- rownames(sum_cont)
sum_cont %>% 
  as.data.frame() %>%
  dplyr::select(varname, mean, sd, median, min, max) %>%
  xtable::xtable(digits = 2) %>%
  xtable::print.xtable(booktabs = TRUE,
                       include.rownames = F)

table(illiberal$party) %>% prop.table()
table(illiberal$gender) %>% prop.table()
table(illiberal$educ) %>% prop.table()
table(illiberal$polint_cat) %>% prop.table()






